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Abstract: It is well known that reflected signals from Global Navigation Satellite Systems 
(GNSS) can be used for altimetry applications, such as monitoring of water levels and 
determining snow height. Due to the interference of these reflected signals and the motion of 
satellites in space, the signal-to-noise ratio (SNR) measured at the receiver slowly oscillates. 
The oscillation rate is proportional to the change in the propagation path difference between 
the direct and reflected signals, which depends on the satellite elevation angle. Assuming 
a known receiver position, it is possible to compute the distance between the antenna and 
the surface of reflection from the measured oscillation rate. This technique is usually 
known as the interference pattern technique (IPT). In this paper, we propose to normalize 
the measurements in order to derive an alternative model of the SNR variations. From 
this model, we define a maximum likelihood estimate of the antenna height that reduces 
the estimation time to a fraction of one period of the SNR variation. We also derive 
the Cramer-Rao lower bound for the IPT and use it to assess the sensitivity of different 
parameters to the estimation of the antenna height. Finally, we propose an experimental 
framework, and we use it to assess our approach with real GPS LI C/A signals. 
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1. Introduction 

Global Navigation Satellite Systems reflectometry (GNSS-R) is a well-established method for 
remotely sensing many relevant geophysical properties of the reflection surfaces. GNSS-R was first 
proposed within the frame of the PAssive Reflectometry Interferometric System (PARIS) project as a 
bistatic radar remote sensing technique for ocean altimetry using the L-band GPS signal [1]. Since then, 
GNSS-R has been demonstrated to be useful in other applications, such as monitoring water levels and 
snow height with a ground approach [2-4]. In this approach, the antenna, situated on a mast, receives 
a direct GNSS signal coming from the satellite and a nadir signal reflected by the observed surface. 
Assuming that the antenna position is known, we can compute the position of the surface of reflection. 
This approach provides precise localization and dating of the measures that allows for spatio-temporal 
comparison of water levels and snow cover, respectively [5-8]. These parameters are very important 
for flood monitoring and avalanche prevention, as well as for hydroelectric companies. Furthermore, 
the approach is noninvasive and can be easily implemented on a portable instrument and embedded in a 
vehicle with a mast [9] . 

GNSS-R altimetry can be carried out in two different ways, depending on the ranging principle: 
code altimetry and phase altimetry [10]. With code altimetry, only the GNSS code delay difference 
between the direct and reflected signals is used. With phase altimetry, the phase of the signal is also 
used for computing this delay difference [9]. The interference pattern technique (IPT) considers the 
behavior of the signal-to-noise ratio (SNR) of the received GNSS signal as a function of the satellite 
elevation [2,11]. The direct and reflected GNSS signals are combined at the antenna. Due to their 
different phase variations, the SNR oscillates at a rate proportional to the distance between the antenna 
and the surface of specular reflection. Unlike satellite or airborne reflection scenarios, ground GNSS 
receivers observe a coherent interference pattern if we consider a flat reflecting surface (compared to 
the carrier wavelength) on an area corresponding to the first Fresnel zone. A few previous works 
can be found analyzing the accuracy of these GNSS-R altimetry techniques [12-15]. Initial works 
proposed simple analytical models to describe the altimetry precision as a function of system/instrument 
parameters [12]. These methods rely on a considerable number of assumptions that might hold only for 
some specific scenarios and provide a pessimistic bound on the achievable precision. In [13], the authors 
proposed a Cramer-Rao lower bound (CRLB) closed expression for code altimetry. The CRLB is a 
statistical tool that provides a lower bound on the achievable estimation error for any unbiased estimator. 
A new derivation was proposed in [15] to compute the CRLB for code altimetry and a specific set of 
measurement data under multiple SNR scenarios. 

Unfortunately, one of the main drawbacks of the IPT is that very long measurement times are usually 
needed. The observed SNR oscillates with the variation of the satellite elevation, but satellite elevation 
varies slowly, thus long measurement times are required to estimate the SNR frequency of oscillation. To 
reduce the estimation time to a fraction of one period of the SNR variations, we propose an alternative 
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model for the measured SNR. This normalized model is based on the normalization of the measured 
signal amplitudes and is possible only after an initial calibration step. This calibration step consists of 
varying the antenna height of the receiver a value dh in order to obtain the minimum and maximum value 
of SNR for a given satellite elevation. Using the normalized model, we define a maximum likelihood 
estimate of the antenna height that allows for the reduction of the required estimation time to a fraction 
of one period of the SNR variation. We also derive the minimum antenna variation range dh as a function 
of the satellite elevation and deduce from this function the minimum observation time as a function of 
the satellite elevation rate. In addition, we derive the CRLB for the IPT and use it to assess the sensitivity 
of different parameters to the estimation of the antenna height. Finally, we propose a novel experimental 
framework, which we use to assess our approach with real signals. 

This paper is organized as follows: Section 2 describes the interference pattern problem. The 
considered signal model is introduced as a function of the receiver height. The proposed estimator is 
presented in Section 3. Section 4 includes the derivation of the CRLB for the IPT and the proposed 
estimator performance assessment using synthetic signals. In Section 5, the proposed experimental 
framework is described, and the results obtained with real GPS LI C/A signals within this framework 
are presented. Finally, Section 6 summarizes the paper, highlighting its main conclusions. 

2. Interference Pattern Problem 

We show in Figure 1 the reflectometry principle for an antenna situated on a mast of height h. In our 
approach, we take into account the subset of satellites M that have a specular reflection on the surface 
to analyze. The antenna receives M scaled, time-delayed and Doppler- shifted signals with known signal 
structures. Each signal corresponds to the line-of-sight or direct signal (sdJ plus its corresponding 
specular reflection (s Ri ). The overall received signal can be modeled as: 

M 

i=i 

M 

= ^2 A DMt - tdJ cos {2ir{f RF + f di )t + (j) D .) 
i=i 

M 

+ ^AR,Ci{t - tr.) cos (27t(/hf + + 4> Di + M)) +n(t), (1) 
i=i 

where f RF is the carrier frequency, fa the Doppler frequency shift of the i-th satellite, <p Di the receiver 
clock phase offset, 4> R . (t) the phase delay between the direct and reflected signals as a function of time, 
TDi,T Ri the time-delays, the pseudorandom code sequence, A F>i and A Ri the amplitudes of the 
direct and reflected received signals and n(t) zero-mean additive Gaussian noise with variance a\. In 
this paper, we will assume that the variations of A F)i and A R . are negligible during the short periods 
of observation considered, e.g., a few minutes. For long observation times, Ar >i and A Ri will change 
with time as a function of the satellite elevation due to the variation of the received power and the 
antenna footprint. The time dependence of td { , r Ri , fa and 0^ has been neglected for simplicity in 
Equation (1), since their variation over time will be compensated for by the receiver's tracking stage 
with little impact on the proposed analysis. 
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Figure 1. Reflectometry principle for an antenna on a mast. 




In our current study, we will consider only the processing of the GPS LI C/A signals. In this case, 
for an antenna mounted on a mast a few meters above the reflecting surface, the difference in GNSS 
signal path between the direct and the reflected signals will be small compared to the chip duration of 
the code. Thus, we can assume that Ri(t — r Ri ) ps Ri(t — r Di ), where Ri(r) is the autocorrelation of 
Cj(i). According to the geometry depicted in Figure 1, it is easy to show that: 

4:71 

Mt) = Y- h M^eiM) (2) 

Ail 

tan(0 di (*)) 

where 9 e ii(t) is the elevation of the i-th satellite at instant t, Xli = 19.042 cm is the GPS LI C/A carrier 
wavelength, h is the height of the antenna and c represents the speed of light. In the following, for 
notation simplicity, we will drop the satellite index i, since we will focus on the processing of the direct 
and reflected signals coming from a single satellite. In this case, by using some trigonometry, we can 
express x(t) as: 



x(t) = A G (t)c (t - r) cos (2tt {{f RF + f d )t) - <f> G (t)) (4) 

where: 

<Pg it) = <f) D + arctan - — — — (5) 

V^ D + A R cos(<j) R (t)) J 



A G (t) = yJA 2 D + A 2 R + 2A D A R cos(0 fl (t)) (6) 

and Aq{€) represents the magnitude of the composite signal, while </>g(£) represents its phase. In practice, 
in a GNSS receiver, the SNR is estimated after the signal down-conversion and correlation with a local 
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code replica. In this case, the SNR is proportional to A 2 G {t), the squared amplitude of the received 
signal. In this context, we see from Equations (2) and (6) how A G {t) evolves as a cosine of the sine 
of the satellite elevation. The frequency of this cosine, -j^, is proportional to the antenna height. This 
means that by estimating the frequency of the observed SNR, we can obtain the height of the receiver. 

In order to get an accurate estimate of the frequency of cos(0#(£)) with classic approaches, one 
must observe at least one period of the signal. For a given initial elevation 9 e i(t 0 ), we define A9 d , the 
satellite elevation variation required to observe one period of the signal. Based on Equation (2) and using 
trigonometric identities, we can thus write: 



— 2 sin 



A9 P 



COS 6ei(t 0 ) + 



A6 e 



0 



(7) 



Figure 2 shows the corresponding elevation variation required according to antenna height, for 
different 9 e i(t 0 ) values. In particular, we can see that for a height of 3 m, one period of the cosine 
can be observed for a satellite elevation variation of at least 2° when 6 e i(to) is close to 0°. If we consider, 
for example, a mean satellite elevation speed u e [ = 1CT 3 °/s, we must thus wait at least 33 min to 
observe one period of the signal. In the next section, we propose a normalization procedure to decrease 
the required observation period. 

Figure 2. Satellite elevation variation as a function of the antenna height. 
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3. Proposed Approach 



3.1. Normalization of the GNSS Signal Amplitudes 

As described in Equation (2), the phase is a function of the satellite elevation and of the antenna 
height. Since the satellite elevation evolves slowly, we propose a calibration procedure that uses instead 
a variation of the antenna height. From Equation (6), the minimum and maximum values of Ac(t) can 
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be obtained when cos((p R {t)) is equal to —1 or one, respectively. They are defined by the following 
expressions: 



(A 
(A 



G.min ) 



A 2 + Al - 2A D A, 



G.max ) 



D -r n R - ^r^ D I\ R 

l + A 2 R + 2A D A R 



(8) 
(9) 



From these two equations, we can deduce the sum of the square of the amplitudes and their product as: 

A 2 + A 2 

a 2 I a 2 G,max ~ G,min ^1A\ 

D ' R 2 ^ ' 

A 2 - A 2 

0 a a G,max G,min 

^s\ D A R — (llj 



Therefore, upon substituting Equations (11) and (10) into Equation (6), the single unknown parameter 
to estimate will be the phase delay <j) R (t), which is proportional to the height of the antenna and the sine 
of the known satellite elevation angle. 

In order to always get the maximum and minimum value of A G , the variation of (fi R should be greater 
than or equal to 2n. According to Equation (2), the minimum variation of the antenna height dh should 
thus be equal to: 



dh r , 



A 



2sin(0 eI ) 



(12) 



In Figure 3, we show the value of dh rnin as a function of the satellite elevation. From this figure, we 
can see that a variation in the antenna height of 0.5 m is sufficient to observe a maximum and a minimum 
value for A G when the satellite elevation is higher than 12°. 

Figure 3. Antenna variation as a function of the satellite elevation. 
1.1 
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3.2. Frequency Estimation 

After down-conversion, the received signal is correlated with a local code replica. For the following 
derivation, we define the samples y[n] as the noisy post-correlation measurements obtained every 
t = nT int , where T int is the coherent integration time. From Equation (6), we can define y[n] as: 

y[n) = A G [n] + w[n] (13) 



A 2 D +A R + 2A D A R cos(<f> R [n])+w[n] (14) 

where w[n] is the resulting zero-mean additive white Gaussian noise (AWGN) with power a 2 . Note 
that if the pre-correlation noise samples are colored (e.g., due to front-end imperfections or interfering 
signals), we assume that spectral whitening has been used to optimally process the RF front-end output 
samples (see, e.g., [16]). From Equation (2), we note that (f) R [n\ evolves linearly as a function of 
sin(6 e i[n]), with a constant factor f3 = An^-h. 
Let us define: 



fmodel 



n}=/3 sm{e*[n]) (15) 



Therefore, the factor (3 defines the frequency of cos(4> R [n}), and its evolution is defined as a function of 
the sine of the elevation. The satellite elevation 8 e i [n] is obtained from the current GPS ephemeris data 
and the estimated position of the GNSS receiver. In order to estimate (3, after calibration, we can define 
the following model for A G [n] : 



I A 2 + A 2 A 2 — A 2 

a r i / G.max ~ G.min , G.max G.min / lmnrtpH "IN rir\ 

a g M = y 2 1 2 cos \ n \> ( 6 ) 

and derive the maximum likelihood estimate of (3 for N measurements as: 

( N 2] 

f3 = argmin < (y[n] — Aa[n] j > (17) 



,n=l 

Finally h is a function of /3 defined by 

ft=/f (18) 

In the next section, we will derive the CRLB for Equation (14) in order to make a feasibility study 
and assess the expected performance of the proposed approach. 

4. Performance Assessment 

4.1. Cramer-Rao Lower Bound 



We are interested in assessing the maximum theoretical accuracy that can be obtained when estimating 
the receiver height h. Unfortunately, Equation (14) is highly nonlinear, which makes it difficult to directly 
assess the impact of its different parameters over the estimation error. This nonlinearity is due mainly to 
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the cosine function in the expression and is exacerbated by the presence of the root mean square. Instead, 
we propose to compute the CRLB for the signal model under consideration. The CRLB provides a lower 
bound on the variance of any unbiased estimator and, thus, will allow us to assess the performance of 
our estimator [17]. 

The signal model considered for Ac [n] is provided in Equation (14). In order to provide more 
insightful results, we will express the reflected signal amplitude as A R = aA D , where a represents 
the attenuation coefficient due to reflection, assumed to be real and less than or equal to one. In addition, 
we define 7 [n] = y sin (9 e i [n]). Thus, we obtain: 



y[n;g\ = A G [n;£] +w[n] = A D y/i + a 2 + 2a cos (7 jnj h) + w[n] (19) 

where £ = \A D , a, h] is our unknown deterministic parameter vector and w[n] are zero-mean AWGN 
samples with variance o\. The CRLB for £ can be expressed as [17]: 



var 



(€>) > [I" (* 



(20) 



where I (0) is the Fisher information matrix (FIM). A full account of the derivation of the FIM for the 
considered signal model of y[n; £] can be found in Appendix 1. The obtained FIM is: 



i(0 



W-1 

71=0 

N-l 

Ad ( a + cos (7 N h)) 

n=0 
N-l 

-A D a J2 7 H sin (7 [n] h) 

n=0 



N-l 

A D J2 (" + cos (7 [n] h)) 

n=Q 

(a+cos(7[n]h)) 



4E 

n=0 
N — l 

4 7[n](a+cos(7[n]/i)) sin(~f[n]h) 



N-l 

-A D a J2 7 N sin (7 [ n ] h) 

n=0 

N-l 

A4, n V 1 7N("+cos(7[n]fe)) sva{y[n]h) 

-A D Ot 2_, J^n 

n=Q 



N-l r 



~A D a E 



AWE 

n^0 



7[n] sin(7[n]/i) 

A G [n;C] 



(21) 



The CRLB for /i can be obtained by computing [I 1 (0)] 33 , resulting in: 

var C RB (hj > gNR • 5 (/i, a, A, A0 el ) (22) 

2(T 2 

where SNR^ = tfa2~ 1S m e post-correlator SNR when only the direct signal is received. 
For simplicity, the g() function is used to represent the multiple terms of [I -1 (0)] 33 . 
A6 e i = {8 e i [0] , 9 e i [1] , . . . , 9 ei [N — l]} is the satellite elevation span covered by N measurements with 
0 < 9 d [n] < 90°. 

The purpose of the following discussion is to identify the effects of {h, a, A, A6 e i} parameters through 
simulation. In order to highlight the effect of a specific parameter, the CRLB is computed for different 
values of the parameter of interest, while the rest are kept fixed. We have set A = Ali, corresponding to 
the GPS C/A LI wavelength, a sampling period T s — 1 sample/s and a = y/0?f (water surface scenario 
for a typical smooth sea [18]) in all of the considered cases, unless otherwise specified. 
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4.2. Assessing the Initial Elevation on the CRLB 



Figure 4 shows the (Jcrb{^) = yvarcRB [hj for different satellite initial elevations (9 e i [0]) and 
antenna heights of h = 2 m and h — 15 m. Elevation variations of 3° were covered by N = 600 
observations starting from the different 9 d [0] values. A reference SNR B = 18 dB was considered. The 
CRLB is computed from the probability density function (pdf) of the data observations. The estimation 
accuracy, lower-bounded by the CRLB, is related to the dependency of the data pdf on the unknown 
parameters (£). The CRLB is higher when the dependency between the observations and the parameter 
to estimate is weak. 



Figure 4. Cramer-Rao lower bound (CRLB) for different initial satellite elevation angles 
and receiver heights. 
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Remarks: 

• When 8 e i [0] is close to zero, the CRLB is several times higher than for higher 6 e i [0] values and 
almost independent of h. In order to explain this behavior, we have shown in Figure 5Left the A G 
evolution for different h values over A9 e i = 1° for 9 e i[0] < 1°. The dependency between h and the 
observations becomes smaller when small 9 e i [0] values are considered, even for the same number 
of samples. 

• When A G [n; £] (see Equation (19)) is at its maximum or minimum value, for a small period 
of observation, the signal tends to a constant equal to A^y/l + a 2 + 2a or A + a 2 — 2a, 
respectively. In this region, the signal evolution can be assumed to be monotone with almost null 
variation. The dependency between the observations and h decreases with the absolute value of 
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the slope of A G [n; £] during the period of observation. In this case, the CRLB is as high as the 
slope of the signal is low. 

• We show in Figure 5Right that the frequency of the cosine evolution of A G [n; 0] decreases when 
9 ei tends to 90°. For h = 2 m and an observation interval of 600 s, we observe periodic CRLB 
increases with the satellite elevation when Equation (19) includes a minimum or a maximum. This 
situation corresponds to the peaks that appear in Figure 4 for h = 2 m and 9 e i[0) > 50°. For 
h = 15 m and an observation interval of 600 s, the minimum frequency of Equation (19) is too 
high to have a constant signal evolution in a period of observation. Thus, the CRLB for h — 15 m 
does not present the peaks as when h = 2 m. 

From the observed behavior of the CRLB, we prefer using the proposed estimator with satellite 
elevations between 10° and 70° and antenna heights greater than 2 m to obtain better estimation accuracy. 
For low elevations, the CRLB is indeed high, independent of the antenna height. For elevations close 
to 90°, the CRLB increases for low antenna heights, because just a small portion of the signal Ac[n; £] 
oscillation, considerably less than a period, is observed. As expected, this effect is exacerbated when the 
uj e i is lower, and a shorter evolution of A G [n; £] is observed for the same measurement period. 

Figure 5. Examples of A G evolution for different receiver heights, in the absence of noise: 
(Left) Aq for different elevation angle variations and receiver heights; (Right) A G for 
elevation angle variations between 60° and 90° and receiver heights of 2 m and 15 m. 



9 ,[0]=0°; A/=600 6 ,[0]=60°; A/=600 




4.3. Assessing the Elevation Rate on the CRLB 

Figure 6 shows the a G B.B(h) computed for different elevation variations A6 e i from 0° to 6° with 
a fixed 9 d [0] = 45°. For receiver heights h = 2 m and h = 15 m, 9 d [0] = 45° was selected, 
since it approximately corresponds to a minimum of the A G [n; £] evolution, which corresponds to the 
worst scenario for the estimation of its frequency of oscillation. These A0 ei values were generated by 
assuming constant u ei up to 10 x 1CT 3 °/s for a fixed measurement period of 600 s and a sampling rate 
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of one sample/s. From the figure, we observe that by increasing A6 e i, the variance decreases quickly 
at first, until an entire period of oscillation of the signal model is covered. This relation can be used 
to set the duration of the measurement time for our estimator, since we are interested in achieving the 
maximum possible accuracy with the minimum number of data observations. Unfortunately, the period 
of oscillation of our model depends on the true value of h and co e i. The measurement time required to 
achieve a certain target accuracy will depend on co e i, and in general, for faster co e i, shorter measurement 
periods will be required to achieve a similar estimation accuracy. 

Figure 6. CRLB for different elevation angle variations and receiver heights. 
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4.4. Assessing the Impact of Amplitude Ratio a on the CRLB 

Figure 7 shows the acRB(h) obtained for different receiver heights (h — 2 m and h = 15 m) and 
reflected/direct amplitude ratios (a) versus SNRj> The coefficient a can be interpreted as the square-root 
magnitude of the polarization-dependent reflection coefficient. For right-hand circular polarization 
(RHCP) of the incident signal and left-hand circular polarization (LHCP) of the reflected signal, a can 
be calculated as a function of 9 e i and the dielectric constant of the scattering surface, assuming a smooth 
surface [18]. Two values for a have been selected. The first value, a w = \/(X7, corresponds to sea 
water, with a typical dielectric constant e w = 73.0 + z57.5 for X Li = 0.19 m and an 9 el = 15° [18]. The 
second value, as = \/0.08, corresponds to fresh snow at —2 °C, with Es = 1-45 + i2.76 x 10" 4 for a 
9 d ~ 10° [19,20]. The A8 e i considered covers a satellite elevation variation of 3° with an 9 d [0] = 15°, 
iV = 600 samples and T s = 1 sample/s. 

The figure shows that the estimation error of h is inversely proportional to a. Some small differences 
are observed between different h values. These differences appear due to the differences on the 
oscillation period covered in each case by the interference pattern signal A G [n; £] used for computing 
the CRLB. 
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Figure 7. o-crbW f° r different receiver heights and reflected/direct amplitude ratios (a) 
versus SNR. 
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4.5. Estimator Performance Evaluation with Synthetic Data 

In order to complement the CRLB analysis, we want to assess the performance of the estimator 
proposed in Section 3 with synthetic data generated using real measurements of the satellites elevation 
(6 e i). These data were generated using the signal model defined in Equation (19) with a constant 
SNRd=18 dB for the direct signal and a sampling period of 1 s. The SNR/) value selected is typically 
achieved at the post-correlation stage and matches the SNR corresponding to a carrier-to-receiver noise 
density C/N 0 = 45 dB-Hz with a front-end bandwidth of 3 MHz, sampling at the Nyquist frequency, 
and a coherent integration time of 1 ms. The power ratio between the reflected and direct signals was set 
to oiyy = 0.7. The satellite elevation measurements correspond to the satellites in view at Calais, France 
(50° 55' 14.093" N, 1° 56' 59.44" E), on 25 September 2013. The root-mean- square error (RMSE) 
was computed for 120 observation periods of 10 min, each period starting one minute later than the 
previous one, with N = 600 observations taken at one sample/s. We assume here that we are in a classic 
bistatic configuration depicted in Figure 1, where the receiver antenna is located at h = 2 m above the 
ground. The parameter (f)R is defined by Equation (2), and the height h is sought with a resolution step of 
1 x 10 -3 m in a bounded search space h space = [0, 5] m. 

In Figure 8a,b, we show the satellites' elevations (9 e i) and their elevation rates (co e i), respectively, as 
a function of time. We present in Figure 8c,d the RMSE of h for the proposed estimator as a function 
of time, computed using the Monte Carlo method with 100 realizations. Figure 8c shows the RMSE 
of the satellites reaching an elevation rate close to zero at some instant during the full observation 
interval. Figure 8d shows the RMSE of the satellites with a high elevation rate, except for the end 
of the observation interval. 
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Figure 8. RMSE obtained for height estimation, with the corresponding elevation and 
elevation rate, for several satellites in view on September 25, 2013. 
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b) Elevation variation rates for satellites in view 




20 40 60 80 100 
Time [min], start :25. 09.201 3 11h52m39s 



d) RMSE : h=2 m, W=600, 
SNR D =18dBanda=V0.7 



120 140 
UTC 



Aft . 




20 40 60 80 100 
Time [min], start :28. 09.201 3 11h52m39s 



120 140 
UTC 



We observe in Figure 8c that the RMSE increases when the satellite elevation rate decreases to less 
than 2 x 10~ 3 °/s. In the case of Satellite 6, the RMSE is higher, because the satellite elevation is close 
to 90°. We observe in Figure 8d that the RMSE of Satellites 3 and 22 increases for high elevations. 
This is due to a low elevation rate and, so, a low frequency of the SNR evolution. These results are in 
accordance with the CRLB study in Section 4.1. 

Let us now compute the RMSE of h for different SNR values and observation interval durations. As 
before, the SNR refers to the signal-to-noise ratio for the direct signal, and the power ratio between the 
direct and reflected signals is kept at = 0.7. In Tables 1-3, we report the RMSE obtained with 1000 
realizations of the noisy Aq signal, for observation intervals of 600 s, 300 s and 150 s, respectively. The 
duration of these observation intervals was selected considering that for a satellite elevation of 35° and 
an elevation rate u e i = 6.8 x 1CT 3 °/s (e.g., Satellite 3 at 12h02 UTC, in Figure 8), a complete period of 
signal is observed after 600 s, a half a period after 300 s and a quarter of a period after 150 s. 

From the results in Tables 1-3, we observe that the proposed estimator is consistent and that the 
RMSE increases when the SNR decreases, which was expected. In these tables, we have defined two 
different sets of satellites. The first set includes Satellites 3, 6, 21 and 22 (in bold in the Tables) with 
a high absolute elevation rate co e i > 6 x 1CT 3 °/s. Satellites 3, 21 and 22 have low initial elevations 
between 20° and 35°, and Satellite 6 has a high elevation, superior to 70°. In the second set, we consider 
the Satellites 7, 16 and 18, with a low u e i. Satellite 7 has a low initial elevation of 18°; Satellite 18 has 
an initial elevation of 44°; and Satellite 16 has a high initial elevation of 75°. 
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Table 1. RMSE of the estimated height h, in meters, for an observation interval of 600 s. 
Reference height h = 2m and a.w = \/0?7 



S atellite^^^^^ 


18 dB 


13 dB 


8dB 


Sat 3 at llh57UTC 


0.001 


0.001 


0.027 


Sat 6 at 12h32 UTC 


0.005 


0.016 


0.090 


Sat 21 at 13h52 UTC 


0.001 


0.001 


0.014 


Sat 22 at 12h42 UTC 


0.001 


0.001 


0.045 


Sat 7 at 13h22 UTC 


0.072 


0.131 


0.432 


Sat 16 at 12hl2 UTC 


0.107 


0.143 


0.596 


Sat 18 at 13hl2UTC 


0.098 


0.1302 


0.484 



Table 2. RMSE of the estimated height h, in meters, for an observation interval of 300 s. 
Reference height h = 2 m and aw = x/O/f. 



^^^^ SNR 
S atellite ^^^^ 


18 dB 


13 dB 


8dB 


Sat 3 at llh57 UTC 


0.005 


0.027 


0.152 


Sat 6 at 12h32 UTC 


0.059 


0.080 


0.286 


Sat 21 at 13h52 UTC 


0.001 


0.018 


0.224 


Sat 22 at 12h42 UTC 


0.001 


0.011 


0.152 


Sat 7 at 13h22 UTC 


0.183 


0.245 


0.849 


Sat 16 at 12hl2 UTC 


0.232 


0.313 


1.162 


Sat 18 at 13M2UTC 


0.457 


0.625 


1.558 



Table 3. RMSE of the estimated height h, in meters, for an observation interval of 150 s. 
Reference height h = 2m and a w = \/o7f. 



^^^^ SNR 
S atellite ^^^^^ 


18 dB 


13 dB 


8 dB 


Sat 3 at Hh57 UTC 


0.116 


0.153 


0.681 


Sat 6 at 12h32 UTC 


0.146 


0.198 


0.708 


Sat 21 at 13h52 UTC 


0.067 


0.111 


0.431 


Sat 22 at 12h42 UTC 


0.123 


0.168 


0.644 


Sat 7 at 13h22 UTC 


0.412 


0.553 


1.473 


Sat 16 at 12M2 UTC 


0.685 


0.937 


1.889 


Sat 18 at 13M2UTC 


1.340 


1.570 


2.464 



For Satellites 3, 21 and 22, the RMSE values presented in Tables 1 and 2 approach the CRLB 
values computed in Section 4.1. In these cases, we reach centimeter accuracy due to the high elevation 
rate of these satellites. However, for a period of observation of 300 s and a C/N 0 = 35 dB-Hz, 
we reach just decimeter accuracy, showing the limits of this technique. In Table 3, we see that 
centimeter accuracy is not reached for an observation interval of 150 s, even with C/N 0 = 50 dB-Hz. 
For 150 s, the SNR evolution covers only a short part of the oscillation period during the interval of 
observation and can be considered monotonic. In this context, we observe in Table 3 that the RMSE 
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strongly depends on the considered part of the SNR cosine evolution rather than on the satellite's tu e i 
(e.g., Satellites 3, 6, 21 and 22). 

For Satellite 6, the results shown in Tables 1 and 2 are less accurate despite the fact that its cu e i is 
similar to those of Satellites 3, 21 and 22. This can be explained by the higher elevation of Satellite 
6. This result is indeed in accordance with the CRLB study (Section 4.1), because in this case, the 
SNR evolves with a lower frequency, so we do not cover a full period of the SNR variation over the 
observation interval. 

For Satellites 7, 16 and 18, the results shown in Tables 1-3 are also in accordance with the expected 
accuracy from the CRLB study. In this case, we reach just meter accuracy, because we observe only a 
small fraction of the SNR variation period due to the the low Q d of the satellites. 

5. Experimental Framework 

5.1. Experimental Results 

In order to assess the proposed method, we constructed the following experimental setup to measure 
the height difference between two antennas, as depicted in Figure 9. The height difference is estimated 
with the interferometric approach described in Section 3. The advantage of the proposed assessment 
method is that we can have centimeter knowledge of the system geometry. Next, we derive the 
link between the height difference of the two antennas and the frequency of variation of the GNSS 
signal power. 

Figure 9. Geometry of the experiment. 




Ground 
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For a height difference h between the two antennas, we have (see Appendix 2): 

a = Axcos(9 e i) cos(A9 Az) — hsm(9 e i), (23) 
h = Axtan(^), (24) 

where a is the path difference between the GNSS signals for the two antennas, Ax the distance between 
the two antennas in the ground plane, 9 e i the elevation for the satellite in view and ip the angle between 
the two antennas and the ground plane. The paraxial approximation is assumed, so that the satellite 
signal arrives at both antennas with the same elevation angle, 9 e [. A9a- z is the angle between the vertical 
plane containing the two antennas and the vertical plane containing the satellite and any of the antennas. 
Then, it follows that: 

a cos 9 d cosA9 Az . a , oc , 

- = sm9 el (25) 

n tamp 



cos 2 A6U 2 /. „ f cosA9 Az \\ 

5 h 1 sin 6 ei + arctan 1 (26) 

tan z (p \ \ tan cp J J 

Ksm{9 el + K 0 ) (27) 



with: 



/cos 2 A9 Az 

K = \ - 2 Az +1, (28) 
tan^ (p 

+ / cosA^4, \ 
K 0 = arctan* — - 1 (29) 



tan <p 

where arctan* () is the quadrant- specific inverse of the tangent. Finally, we can write: 

a = h K sin (9 e i + K 0 ) (30) 
The phase delay between the direct signal received by the two antennas, 4>r P , can be expressed as: 

2tt 2tt 2tt f cos 9 ei cosA# A2 \ 2tt . 

4>r = -r-a = — JLia = — f L ih sm9 e i = — f L ih K sin 0 ei + K 0 ) (31) 

Ac c \ tamp J c 

(j)ft P depends on 9 e \ and on the satellite azimuth, with A9a z assumed to be constant over the entire 
observation time. We then can conclude that 0^f p evolves linearly, with a slope (3 exp = —JliK as a 
function of K sin(# e z + K 0 ). Finally, f3 exp = |/3, where 0 is the frequency of the variation of the SNR 
as a function of sin(9 e i) in the reflection scenario described in Section 2. 

5.2. Assessment with Real Data 



Figure 10 shows the experimental vehicle and the telescopic mast used. The figure also shows the 
system that gives us the ability to precisely change the height of the antenna installed on the top of 
the mast. The second direct antenna is situated on the roof of the car at a horizontal distance of 
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Ax = 1.92 m. The height h between the two antennas is known, so we can derive the value of ip 
and define the complete geometry of the experimental system. We used the Novatel OEM4-G2 ProPak 
RT2W (GPS + WAAS/EGNOS) commercial receiver [21]. The two antennas were connected using a 
passive RF-combiner. The C/N 0 measurements, as well as the satellite elevation values were provided 
by the receiver at a rate of 1 Hz. 



We show in Figure 11 an example of C/N 0 evolution as a function of K sm(9 e i + K 0 ). In this 
figure, we differentiate two periods of time in the signal. These two periods correspond to two different 
processing steps: the calibration step and the estimation step. 

We report in Tables 4 and 5 the estimated height (h) obtained at Calais, France (50° 55' 14.093" N, 
1° 56' 59.44" E), on 17 January 2014, with the proposed method. The reference height used was 
href = 2.13 m at 14h09 UTC and h ref = 8.24 m at 14h50 UTC. These heights were manually tape 
measured. Figures 12 and 13 show the constellation of the visible satellites during the measurement 
periods and the direction of the experimental setup (working direction). We plotted in these figures the 
satellites' trajectories, with a star marking the end of each trajectory. In the experiments, h was estimated 
with N = 600 observation samples and an search step resolution of 1 x 10 -3 m in a bounded interval of 
h = [h re f — 2, h re f + 2] m. 

For h re f = 2.13 m, the signals from Satellites 5, 9, 20, 16 and 29 were not considered due to their 
low C/Nq and low elevation angles. In these cases, the proposed estimator performed poorly, and the 
estimated height was far from h re f- The estimated height obtained for Satellites 4, 7, 10 and 23 showed 
a difference with the h re f between 2 cm and 9 cm. The mean estimated height is 2.15 m. 



Figure 10. The vehicle and its telescopic mast. 





Sensors 2014, 14 



10251 



Figure 11. Evolution of C/N 0 for h re j = 2.13 m and Ax = 1.92 m. 

C/N Q evolution of satellite 4 the 1 7 January 201 4 at 1 4h09m44s UTC 



-. 

•: 
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Table 4. 17 January 2014, at 14h09 UTC. Reference height: h = 2.13 m, Ax = 1.92 m. 



Satellite 
PRN 


Estimated 
height (h) (m) 


\K sin(8 el + K 0 )\ 
mean variation (s _1 ) 


(dB-Hz) 


C/No^ max 
(dB-Hz) 


Comment 


2** 


2.86 


7.0 x 10~ 5 


43.8 


50.5 


low variation of \K sin(6* e ; + K 0 )\ 


4 


2.22 


1.7 x 10~ 4 


44.8 


50.2 




7 


2.16 


1.5 x 10~ 4 


45.6 


49.8 




8* 


1.93 


1.6 x 10~ 4 


39.5 


47 


low C/Nq 


10 


2.10 


1.0 x 10~ 4 


47.5 


51.4 




13** 


1.53 


9.6 x 10~ 5 


48.4 


51.2 


low variation of \K sin(# e ; + Kq)\ 


16 


2.11 


1.4 x 10~ 4 


35.5 


45.3 


low C/N 0 


23 


2.19 


1.5 x 10~ 4 


44.1 


49.8 





Table 5. 17 January 2014, at 14h09 UTC. Reference height: h = 8.24 m, Ax = 1.92 m. 



Satellite 
PRN 


Estimated 
height (h) (m) 


\Ksm(9 el + K 0 )\ 
mean variation (s _1 ) 


C/N 0min 
(dB-Hz) 


(dB-Hz) 


Comment 


2 


8.29 


7.2 x 10~ 5 


44.5 


50.2 


low variation of \K sm{6 e i + Ko)\ 


5 


8.27 


1 x 10" 4 


42.1 


49.2 




7** 


8.79 


5.9 x 10~ 5 


47.6 


50.9 


low variation of \K sin(# e ; + K 0 )\ 


8 


8.23 


1.2 x 10~ 4 


43.9 


47.8 




9 


8.19 


1.3 x 10~ 4 


44 


49.1 




10** 


8.92 


5.5 x 10~ 5 


48.4 


51.2 


low variation of \K sin(# e ; + Kq)\ 


13* 


8.01 


8.9 x 10~ 5 


45.6 


50.8 


low variation of \K sin(6* e ; + Kq)\ 
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Figure 12. Trajectories of the satellites in view during the measurements, on 17 January 
2014, at 14h09 UTC. 
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NORTH 
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Figure 13. Trajectories of the satellites in view during the measurements on 17 January 
2014, at 14h50 UTC. 



SATELLITE SKYPLOT : 17.01.2014, start at 14h50m UTC 

NORTH 




'°° experimental 
SOUTH setup at 180° 

For h re f = 8.24 m, the signals from Satellites 16, 23, 26 and 29 were again not considered due to 
their low C/N 0 s. The results obtained with Satellites 2, 5, 8 and 9 are closer to h re f, with a difference 
between 1 cm and 5 cm. In this case, the mean estimated height is 8.25 m. 

We can conclude that, after the calibration step, the proposed estimator can achieve centimeter 
accuracy under the experimental setup for a period of observation of 600 s. 
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6. Conclusions 

In this article, we used an IPT to estimate the height between an antenna and a ground surface, where a 
GNSS signal has been reflected. We proposed to normalize the SNR measurements in order to construct 
a model of its evolution over time. The proposed estimator is based on two steps: A calibration step 
and an estimation step. The aim of the calibration step is to measure the maximum and minimum values 
of the SNR (or, equivalently, the C/N 0 ) amplitude, in order to model the SNR variations for a bounded 
interval of possible surface heights. 

The maximum likelihood estimate of the antenna height constructed with this nonlinear model was 
assessed in a study of the CRLB of the model. In this study, we showed that the accuracy of the estimation 
can be defined as a function of the satellite elevation, the elevation rate, the C/Nq and the power ratio 
between the direct and reflected signal. 

In order to assess the method, we used synthetic data and verified that the proposed estimator is 
consistent with the number of observation samples and the C/N 0 of the GNSS signal. We also showed 
that, using the proposed calibration step, we can expect centimeter accuracy with half a period of the 
SNR oscillation when we are in the best scenario. These conditions were identified using the CRLB. 

Finally, we proposed an experimental framework that uses two direct signals received in different 
locations. The results using real data from field measurements showed that the proposed estimator can 
provide centimeter accuracy for a period of observation of 10 minutes within this framework. 
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Appendix 1: FIM Computation 

The Fisher information matrix (I (£)) for a deterministic parameter vector £ is defined as [17]: 
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-E 



d 2 In p(x;$ > ) 



(32) 



The observed samples are modeled as y [n; £] = A G [n; 0] + w [n], where A G [n; 0} is the signal 
component and w [n] is zero-mean AWGN with variance a 2 . In general, for a vector signal parameter 
estimated in the presence of AWGN, Equation (32) yields: 



[1(0] 



1 N-1 



ds [n;£] ds [n; £] 



n=0 



(33) 



In Section 4.1, the signal model A G [n; £] considered is: 



A G [n; 0] = Ada/1 + a 2 + 2a cos (7 [n] /i) 



(34) 



where £ = [ A D , a, h ] T is the unknown vector parameter. 

Every element of the FIM, is computed using a combination of the following partial 

derivatives of A G [n; 0] : 



ds [n; £] 



1 



dA D A D 

ds[n;g\ A 2 D (a + cos (7 [n] h)) 

da A G [n; £] 

ds [n; £] A 2 D a^ [n] sin (7 [n] h) 

dh Aq [n;£] 

In this way, the FIM is defined as a symmetric 3x3 matrix as follows 

r J n-i J N-i 



(35) 



[n;Z] 



n=0 



AdJ2 ( a + cos (7 M h)) 

2 



n=0 
N-1 



(a+cos(7[n]h)) 



N-1 

A D J2 (a + cos (7 [n] h)) A% £ 

n=0 n=0 

N-1 N-1 

4 7[n](a+cos(7[n]/i)) sin(7[n]/i) 



N-1 

—Aj^a Y^, 7 [ n ] si n (7 [ n ] h) 

n=0 

N — l 

^2 7M("+cos(7[n]fa)) sin(i[n]h) 



n=0 



-A D a J2 7 N sin (7 [n] h) -A 4 D a 

n=0 n=Q 



AW £ 



r ~y[n] sin(7[n]/i) 1 

A G [n;C] 



n=0 



(36) 



Appendix 2: Experimental Framework Geometry Computation 

Figure Al shows the geometry of the experimental framework and defines two frames, (0,x,y,z) and 
(0,X,Y,Z), for the geometry computation. The plane (0,x,z) contains the car antenna A and the mast 
antenna B. The plane (0,X,Z) contains the mast antenna and the satellite in view. 

The relation between both frames is: 

X = x cos A8az + y sin A9a z 

Y = -x sin A6 Az + V cos A9az ( 37 ) 
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The path difference a is defined with point H, which corresponds to the orthogonal projection of point 
A over BD, where BD is the line from the mast antenna to the satellite in view. We thus have Y# = 0, 

so y H = x H tan 6 el and X H = cos ^ Az ■ 

Xh-x B b = tan 6^, with X B = 0 and Z B = h, implies that Z H — h = X H tan 6^/ and 

Z H = zh = XH ^ + h. 

n n cos A0Az 

Figure Al. Geometry of the experimental framework. 




(AH) and (BH) being orthogonal, AH ■ BH = 0, so: 

(x H - Ax)x H + y 2 H + zh{z h - h) = 0 (38) 

/ 9 . „ tan 2 9 e i \ 9 f h tan 0 e i . \ 

l + tan 2 A^ + — f x\ + — ^-Ax ) x H = 0 (39) 

V cos 2 AO Az J \cosA9 Az J 

1 h tan 9 e i 

—- YXfT XH = Ax XfT (40) 

cos 9 e l COS 2 A9az COS A9az 

= Ax cos 2 9 ei cos 2 A9az — h sin 9 e i cos 9 e i cos AO a z (41) 

(42) 



Finally, 



And: 



Ch = t±X COS 0 e l COS Az ~ n Sin U e l COS V e l COS LAUAz 

xh = (Ax cos 0 e i cos AO Az — h sinOei) COS 0 e l COsA^Az- 
||M|| 2 = x 2 H + y 2 H + {z H - h) 2 

±H_ 

cos 2 0 e i cos 2 AOaz 
= (Ax cos 0 e i cos AO Az — h sin 0 e l) 

= \\BH\\ = I As cos 6^ cosAOaz — h sin 0 e i\ (46) 



(43) 
(44) 
(45) 



a 
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